#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# SConstruct  (Madagascar Script)
#
# Purpose: Modeling of a seismic data cube A(m,h,t) using Kirchhoff modeling.
#
# Site: https://www.geofisicando.com
# 
# Programmer: Rodolfo A. C. Neves (Dirack) 18/09/2019
#
# Email: rodolfo_profissional@hotmail.com
#
# License: GPL-3.0 <https://www.gnu.org/licenses/gpl-3.0.txt>.

__author__="Rodolfo Dirack <rodolfo_profissional@hotmail.com>"
__version__="1.0"
__site__="https://www.geofisicando.com"

Help('''
    Kirchhoff modeling example for a gaussian reflector
    ===================================================

    Author: %s
    Version: %s
    Site: %s
    ''' % (__author__,__version__,__site__))

# Madagascar library
from rsf.proj import *

# Python math library
import math

# Ploting functions
def wiggle(title):
    return '''
    wiggle wheretitle=top yreverse=y transp=y label1=Time unit1=s label2=CMP unit2=km pclip=99.5 title="%s" 
    ''' % title

def grey(title):
    return '''
    grey label1=Time unit1=s label2=CMP unit2=km pclip=99.5 title="%s" 
    ''' % title

def grey3(title):
	return '''
	byte |
	transp plane=23 |
	grey3 flat=n frame1=500 frame3=80 frame2=200
	label1=Time unit1=s 
	label3=Offset unit3=km 
	label2=CMP unit2=km
	title="%s" point1=0.8 point2=0.8 
	''' % title

# Modeling: Gaussian reflector in a velocity linear model
# velocity increases with depth with a 0.5 velocity gradient
Flow('gaussianReflector',None,
     '''
     math d1=0.01 n1=2001 o1=-5 unit1=km label1=Offset
     output="4-3*exp(-(x1-5)^2/9)"
     ''')

for g in range(2):
    gaussianReflector = 'gaussianReflector%d' % g
    Plot(gaussianReflector,'gaussianReflector',
         '''
         graph min2=0 max2=4 min1=0 max1=10
         yreverse=y plotcol=%d plotfat=%d
         wantaxis=n wanttitle=n scalebar=y pad=n
         ''' % ((7,0)[g],(7,3)[g]))

# Velocity Model
Flow('velocityModel','gaussianReflector',
     '''
     window min1=0 max1=10 |
     spray axis=1 n=451 d=0.01 o=0 label=Depth unit=km |
     math output="1.5+0.5*x1+0.0*x2"
     ''')

Plot('velocityModel',
     '''
     grey color=j allpos=y bias=1.5 scalebar=y wanttitle=n
     barreverse=y barlabel=Velocity barunit=km/s
     ''')

Result('gaussianReflectorVelocityModel','velocityModel gaussianReflector0 gaussianReflector1','Overlay')

Flow('reflectorDip','gaussianReflector','math output="2/3*(x1-5)*input" ')

# Kirchhoff Modeling
Flow('dataCube','gaussianReflector reflectorDip',
     '''
     kirmod cmp=y dip=${SOURCES[1]} 
     nh=161 dh=0.025 h0=0
     ns=401 ds=0.025 s0=0
     freq=10 dt=0.004 nt=1001
     vel=1.5 gradz=0.5 gradx=0.0 verb=y |
     put d2=0.0125 label3="CMP" unit3="Km" label2="Offset" unit2="Km" label1=Time unit1=s
     ''')

Result('dataCube',grey3('Modeled seismic data cube'))

# Select a CMP gather m0=5Km
Flow('cmpGather','dataCube',
	'''
	window n3=1 f3=200
	''')

Result('cmpGather',grey('CMP Gather 5Km'))

End()
